This document contains the code tutorial from “A comparison of statistical models used to characterize species-habitat associations with movement data”. This code tutorial demonstrates the analysis of ringed seal movement in eastern Hudson Bay, Canada.
library(tidyverse)
library(ggplot2)
library(lubridate)
library(rnaturalearth)
library(rgdal)
library(terra)
library(amt)
library(momentuHMM)
library(viridis)
library(sf)
library(here)Next we load in the fish data. This is a subset of the data from Florko et al. (2021).
# load fish data
fish <- read.csv(here("data/preydiv.csv"))Visualize the fish data.
#prepare world data for mapping
## List of azimuthal equal area - HudsonBay
natearth <- ne_countries(scale = "medium",returnclass = "sp")
natearth <- natearth[which(natearth$continent!="Antarctica"),]
nat_trans <- spTransform(natearth,CRS("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))## Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
## to PROJ string
## Warning in wkt(obj): CRS object has no comment
# plot fish map
fishmap <- ggplot() +
geom_tile(data = fish, aes(x = lon,y = lat,fill = preydiv)) +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
fishmapNext we load in the seal movement data. This is a subset of the data from Florko et al. (2023).
# load seal data
seal <- read.csv(here("data/seal_track_m.csv"))
head(seal)## lon lat date id
## 1 357940.9 -368949.6 2012-10-29 14:00 116484_1
## 2 371594.3 -353092.9 2012-10-30 14:00 116484_1
## 3 388429.2 -361478.4 2012-10-31 14:00 116484_1
## 4 399326.2 -357288.1 2012-11-01 14:00 116484_1
## 5 411788.2 -349500.7 2012-11-02 14:00 116484_1
## 6 407708.1 -372564.6 2012-11-03 14:00 116484_1
# ensure the data is in the correct format
seal <- seal %>%
mutate(id = as.character(id),
date = as.Date(date)) Visualize the seal data on top of the fish data.
# plot p
sealfishmap <- fishmap +
geom_point(data=seal, aes(x=lon, y=lat), alpha = 0.6, color = "#FCEEAE") +
geom_path(data=seal, aes(x=lon, y=lat), color = "#FCEEAE")
sealfishmapWe will fit four models: 1) a resource selection function (RSF), 2) a step selection function (SSF), 3) a SSF that allows the habitat covariate to affect movement, and 4) a hidden Markov model (HMM). All four of these models will include prey diversity as a covariate.
We will fit the RSF first, as it is the most simple of the four models. In preparation, we will generate an availability sample, create a map to visualize the spatial distribution of the used and available locations, and extract covariates for the used and available locations. Next, we will fit the model and view the model summary.
All of the step selection functions (both the RSF and two SSFs) will be fit using the ‘amt’ package (Signer et al. 2019).
# prep data and generate availability sample
set.seed(2023)
data_rsf <- seal %>%
make_track(lon, lat, date) %>% # convert data to track format
random_points() # generate availability sample; default is ten times as many available points as observed points
# plot used vs available locations on-top of prey diversity
sealfishmap +
geom_point(data=data_rsf, aes(x=x_, y=y_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type") # rasterize and extract prey diversity covariate
fish_raster <- terra::rast(fish, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
data_rsf <- data_rsf %>%
extract_covariates(fish_raster)
# fit rsf (a binomial logistic regression)
rsf1 <- data_rsf %>%
amt::fit_rsf(case_ ~ preydiv, model = TRUE)
# see model summary
summary(rsf1)##
## Call:
## stats::glm(formula = formula, family = stats::binomial(link = "logit"),
## data = data, model = TRUE)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5753 -0.4729 -0.4294 -0.3754 2.3363
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.118 1.400 -4.369 1.25e-05 ***
## preydiv 6.353 2.312 2.748 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 938.28 on 1539 degrees of freedom
## Residual deviance: 930.60 on 1538 degrees of freedom
## AIC: 934.6
##
## Number of Fisher Scoring iterations: 5
Next we will fit our two SSFs. The workflow is similar to that of the RSF.
# prep data and generate availability sample
set.seed(2023)
data_ssf <- seal %>%
make_track(lon, lat, date) %>% # convert data to track format
steps() %>% # convert track data to step format (i.e., with a start and an end)
random_steps() %>% # generate availability sample
arrange(case_)
# plot used vs available locations on-top of prey diversity
sealfishmap +
geom_point(data=data_ssf, aes(x=x2_, y=y2_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type") # extract prey diversity covariate
data_ssf <- data_ssf %>%
extract_covariates(fish_raster, where = "end") #sample at end of step
# transform movement covariates
data_ssf <- data_ssf %>%
mutate(cos_ta = cos(ta_),
log_sl = log(sl_))
# fit ssfs
## ssf1: ssf without covariate affecting movement kernel
ssf1 <- data_ssf %>%
fit_clogit(case_ ~ preydiv + log_sl + cos_ta + strata(step_id_), model = TRUE)
## ssf2: ssf with covariate affecting movement kernel
ssf2 <- data_ssf %>%
fit_clogit(case_ ~ preydiv*log_sl + cos_ta + strata(step_id_), model = TRUE)
# see model summaries
summary(ssf1)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv + log_sl +
## cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## preydiv 0.957799 2.605953 6.772711 0.141 0.888
## log_sl -0.008404 0.991631 0.083328 -0.101 0.920
## cos_ta 0.031707 1.032215 0.130643 0.243 0.808
##
## exp(coef) exp(-coef) lower .95 upper .95
## preydiv 2.6060 0.3837 4.477e-06 1.517e+06
## log_sl 0.9916 1.0084 8.422e-01 1.168e+00
## cos_ta 1.0322 0.9688 7.990e-01 1.333e+00
##
## Concordance= 0.494 (se = 0.029 )
## Likelihood ratio test= 0.09 on 3 df, p=1
## Wald test = 0.09 on 3 df, p=1
## Score (logrank) test = 0.09 on 3 df, p=1
summary(ssf2)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv * log_sl +
## cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## preydiv 2.791e+01 1.318e+12 2.177e+01 1.282 0.200
## log_sl 1.663e+00 5.273e+00 1.286e+00 1.293 0.196
## cos_ta 2.338e-02 1.024e+00 1.311e-01 0.178 0.858
## preydiv:log_sl -2.744e+00 6.433e-02 2.101e+00 -1.306 0.192
##
## exp(coef) exp(-coef) lower .95 upper .95
## preydiv 1.318e+12 7.589e-13 3.919e-07 4.430e+30
## log_sl 5.273e+00 1.897e-01 4.241e-01 6.556e+01
## cos_ta 1.024e+00 9.769e-01 7.917e-01 1.324e+00
## preydiv:log_sl 6.433e-02 1.555e+01 1.046e-03 3.955e+00
##
## Concordance= 0.555 (se = 0.029 )
## Likelihood ratio test= 1.81 on 4 df, p=0.8
## Wald test = 1.81 on 4 df, p=0.8
## Score (logrank) test = 1.81 on 4 df, p=0.8
Finally, we will fit the HMM using momentuHMM (McClintock and Michelot, 2018). In preparation, we define initial parameters, and then update the parameters using our fit model, to ultimately fit a more refined model.
# prepare dataset
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame()
# convert into momentuHMM format
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("x", "y"),
type = "UTM",
covNames = c("preydiv"))
# define initial parameters
nbStates <- 3 # number of states
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5000, 18000, 38000) # mean step length for each state
sigma0 <- c(5000, 10000, 10000) # sd step length for each state
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5) # turning angle for each state
formula = ~ preydiv # identify covariates
# fit basic 3-state hmm
set.seed(2023)
hmm1 <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0))
# retrieve parameters to refine model
Par0_hmm1 <- momentuHMM::getPar0(hmm1, formula=formula)
# fit a refined hmm with parameters from hhm1
set.seed(2023)
hmm2 <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm1$Par,
delta0 = Par0_hmm1$delta,
beta0 = Par0_hmm1$beta,
formula=formula)We plot the decoded states (estimated behaviours).
# plot decoded states
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2))
hmmstate_plot <- ggplot() +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_path(data=data_hmm, aes(x=x, y=y, color = state, group =ID)) +
geom_point(data=data_hmm, aes(x=x, y=y, color = state, shape = state), size=2, alpha = 0.8) +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"),
labels = c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"),
name = "HMM state") +
scale_shape_manual(values = c(15, 16, 17),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"),
name = "HMM state") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
hmmstate_plotPlot a histogram of step lengths for each state. In order to do this, we need to transform the movement data to UTM and refit our HMM in order to get step length in a meaningful unit (i.e., kilometers)
## prep data in lat/lon and refit model (so step length is in kilometers)
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y")) %>%
st_set_crs("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") %>%
st_transform("+proj=longlat +datum=WGS84") %>%
mutate(long = unlist(map(geometry,1)),
lat = unlist(map(geometry,2))) %>%
st_drop_geometry()
# do momentuhmm which calculates step length in KM
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("long", "lat"), type = "LL", covNames = c("preydiv"))
# define parameters
nbStates <- 3
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5, 12, 38)
sigma0 <- c(3, 5, 8)
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5)
# fit HMM (with step in km)
set.seed(2023)
hmm1_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0))
# get parameters
Par0_hmm1_km <- momentuHMM::getPar0(hmm1_km, formula=formula)
# fit HMM with parameters from hmm1_km
set.seed(2023)
hmm2_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm1_km$Par,
delta0 = Par0_hmm1_km$delta,
beta0 = Par0_hmm1_km$beta,
formula=formula)
# add the state estimate from the HMM
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2_km))
# calculate frequencies of states
v <- momentuHMM::viterbi(hmm2_km)
stateFreq <- table(v) / length(v)
# plot colours
colours.states <- c("#99DDB6", "#539D9C", "#312C66")
# generate sequence for x axis of density functions
x <- seq(0, 50, length=1000)
# get converged mean and sd for each state
meanARS <- hmm2_km$mle$step[1,1]
sdARS <- hmm2_km$mle$step[2,1]
meanCR <- hmm2_km$mle$step[1,2]
sdCR <- hmm2_km$mle$step[2,2]
meanTR <- hmm2_km$mle$step[1,3]
sdTR <- hmm2_km$mle$step[2,3]
# calculate shape and scale of the gamma distributions from mean and sd
sh <- function(mean, sd) { return(mean^2 / sd^2)}
sc <- function(mean, sd) { return(sd^2 / mean)}
# get density functions of the distributions
y_ARS <- dgamma(x, shape=sh(meanARS,sdARS), scale=sc(meanARS,sdARS)) * stateFreq[[1]]
y_CR <- dgamma(x, shape=sh(meanCR,sdCR), scale=sc(meanCR,sdCR)) * stateFreq[[2]]
y_TR <- dgamma(x, shape=sh(meanTR,sdTR), scale=sc(meanTR,sdTR)) * stateFreq[[3]]
# combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging", x=x)
df.y_CR <- data.frame(dens=y_CR, State="ARS", x=x)
df.y_TR <- data.frame(dens=y_TR, State="Travelling", x=x)
statedis <- rbind(df.y_ARS, df.y_CR, df.y_TR)
# plot distributions
hmmstepdist_plot <- ggplot() +
geom_line(data=statedis,aes(x=x,y=dens,colour=State,linetype=State), size=1.2) +
scale_colour_manual(values=c(colours.states,"#000000"),
breaks = c('Foraging', 'ARS', 'Travelling'),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
scale_linetype_manual(values=c("solid","solid", "solid"),
breaks = c('Foraging', 'ARS', 'Travelling'),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
ylab("Density") +
xlab("Step length (kms)") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## Warning: Please use `linewidth` instead.
hmmstepdist_plotPlot a histogram of turning angle for each state.
# generate sequence for x axis of density functions
x <- seq(-pi, pi,length=1000)
# get converged mean and concentration for each state
meanARS <- hmm2_km$mle$angle[1,1]
sdARS <- hmm2_km$mle$angle[2,1]
meanCR <- hmm2_km$mle$angle[1,2]
sdCR <- hmm2_km$mle$angle[2,2]
meanTR <- hmm2_km$mle$angle[1,3]
sdTR <- hmm2_km$mle$angle[2,3]
# get density functions of the distributions
y_ARS <- CircStats::dvm(x, mu=meanARS, kappa=sdARS) * stateFreq[[1]]
y_CR <- CircStats::dvm(x, mu=meanCR, kappa=sdCR) * stateFreq[[2]]
y_TR <- CircStats::dvm(x, mu=meanTR, kappa=sdTR) * stateFreq[[3]]
# combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging",x=x)
df.y_CR <- data.frame(dens=y_CR, State="ARS",x=x)
df.y_TR <- data.frame(dens=y_TR, State="Travelling",x=x)
cmb <- rbind(df.y_TR, df.y_CR, df.y_ARS)
# plot distributions
hmmangledist_plot <- ggplot() +
geom_line(data=cmb,aes(x=x,y=dens,colour=State), size = 1.2) +
scale_colour_manual(values=c(colours.states), breaks = c('Foraging', 'ARS', 'Travelling'), labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
scale_x_continuous(limits=c(-pi,pi))+
ylab("Density") +
xlab("Turn angle (radians) ") +
theme_minimal()
hmmangledist_plotNow we will estimate the utilization distribution from each model to demonstrate the inference from each model. The utilization distribution is defined as the two-dimensional relative frequency distribution of space use of an animal (Van Winkle 1975). This is a simple calculation for the RSF, where we multiply the model coefficent with the resource (prey diversity), exponentiate (since it is a logistic regression), and normalize the estimate. The calculations are more complex for the SSFs since they are conditional models that integrates the movement process. Thus, for the SSFs we calculate the steady-state utilization distribution (SSUD), which is the long-term expectation of the space-use distribution across the landscape (Signer et al. 2017). Luckily, amt has functions to estimate the SSUD.
First, we prepare a dataframe to predict on.
# prep the fish data
newfish <- fish_raster %>%
terra::as.data.frame(xy = TRUE) %>%
filter(x > 100000 & x < 600000 & y > -550000 & y < 0)Now we will predict the estimated probability of use from the RSF by hand.
# grab model coefficients
modcoef <- summary(rsf1)$coef
# prediction for each cell
x <- exp(modcoef[1] + (modcoef[2] * newfish$preydiv))
# scale the data
x <- x / (1 + x)
# range fn
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
# set the range from zero to one
newfish$rsf_prediction <-range01(x)
# plot
map_RSF <- ggplot() +
geom_tile(data = newfish, aes(x = x,y = y, fill = rsf_prediction)) +
scale_fill_viridis(option = "mako", name = "RSF prediction", limits = c(0, 0.3)) +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygons
map_RSFWe will first estimate the SSUDs from the simple SSF that does not allow prey diversity to affect the movement kernel. We need to refit our SSF model to include terms for the home range in order to make sure our simulated track stays within our study area.
# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date) %>%
steps() %>%
random_steps() %>%
arrange(case_) %>%
amt::extract_covariates(fish_raster, where = "both") %>%
na.omit()
# fit SSF1 model
m1 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + cos(ta_) + log(sl_) +
# add terms for home range
x2_ + y2_ + I(x2_^2 + y2_^2) +
strata(step_id_))We will now simulate a track and visually observe it.
# set starting position for the simulation
set.seed(2023)
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))
# random sample along the track as the start
[sample(1:nrow(seal), 1), ],
dt_ = hours(2))
# generate redistribution kernel
k1 <- redistribution_kernel(m1, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 0.1,
n.control = 1e3)
# Now simulate a path of length 1e3
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p1 <- amt::simulate_path(x = k1, n = n_steps, start = start, verbose = TRUE)
# plot simulated track
ssf_track_1 <- fishmap +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
geom_point(data = p1, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p1, aes(x = x_,y = y_))## Regions defined for each Polygons
ssf_track_1We can see that it mostly stays within the study area. We will use this track to estimate the SSUD, and visualize the results.
# estimate SSUD
uds_ssf1 <- tibble(rep = 1:n_steps1,
x_ = p1$x_, y_ = p1$y_,
t_ = p1$t_, dt = p1$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "global") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)
# plot SSUD
map_SSF1 <- ggplot() +
geom_tile(data = uds_ssf1, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF1 prediction") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygons
map_SSF1We will follow the same steps to generate a SSUD for the SSF that allows prey diversity to affect the movement kernel.
# fit SSF2 model
m2 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + cos(ta_)+
preydiv_end:log(sl_) +
# add terms for home range
x2_ + y2_ + I(x2_^2 + y2_^2) +
strata(step_id_))
# set starting position for the simulation
set.seed(2023)
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))
# first location of the track as the start
[sample(1:nrow(seal), 1), ],
dt_ = hours(2))
# generate redistribution kernel
k2 <- redistribution_kernel(m2, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 0.3,
n.control = 1e3)
# Now simulate a path of length 1e3
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p2 <- amt::simulate_path(x = k2, n = n_steps, start = start)
# plot simulated track
ssf_track_2 <- fishmap +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
geom_point(data = p2, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p2, aes(x = x_,y = y_)) ## Regions defined for each Polygons
ssf_track_2# estimate SSUD
uds_ssf2 <- tibble(rep = 1:n_steps1,
x_ = p1$x_, y_ = p1$y_,
t_ = p1$t_, dt = p1$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "local") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)
# plot SSUD
map_SSF2 <- ggplot() +
geom_tile(data = uds_ssf2, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF2 prediction") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygons
map_SSF2We will first estimate the stationary state probabilities of each state based on prey diversity. This is easily done using momentuHMM’s function stationary().
# grab estimated stationary state probabilities from our fitted HMM
x <- as.data.frame(momentuHMM::stationary(hmm2, data.frame(preydiv = newfish$preydiv)))
newfish$hmm_state1 <- x$state.1
newfish$hmm_state2 <- x$state.2
newfish$hmm_state3 <- x$state.3
# prepare data
newfish_long <- newfish %>%
tidyr::pivot_longer(cols = hmm_state1:hmm_state3,
names_to = "model", values_to = "prediction") %>%
mutate(dplyr::across(model, factor, levels=
c("hmm_state1", "hmm_state2", "hmm_state3")))
# plot
map_HMM <- ggplot() +
geom_tile(data = newfish_long, aes(x = x, y = y, fill = prediction)) +
scale_fill_viridis(option = "mako", limits = c(0,1)) +
labs(fill = 'HMM predicted\nprobability') +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), colour = "white", fill = "grey90", size = 0.3) +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
facet_wrap(~model, labeller = as_labeller(c('hmm_state1' = "Localized search",
'hmm_state2' = "Area-restricted search",
'hmm_state3' = "Travel")))## Regions defined for each Polygons
map_HMMbase <- newfish %>%
# (note, it doesn't matter what values you pick for sl_ and ta_)
mutate(log_sl = log(45),
cos_ta = cos(1))
x1 <- base
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(rsf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
# PLOT
rsf_line <- ggplot(res, aes(x = preydiv_x1, y = (log_rss))) +
geom_line(size = 1, color = "#F8D59F",linetype = 2) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.4, fill = "#F8D59F") +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()## Prediction for SSF1
x1 <- base
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
# PLOT
ssf_line_1 <- ggplot() +
geom_line(data = res, aes(x = preydiv_x1, y = log_rss), size = 1, linetype = 3) +
geom_ribbon(data = res, aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.3) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
## Prediction for SSF2
nums = seal %>%
make_track(lon, lat, date) %>%
steps %>%
summarize(quants = quantile(sl_, c(0.25, 0.5, 0.75))) %>%
pull()
x1 <- base
results <- lapply(nums, function(i) {
x1$log_sl <- i
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf2, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
} )
results <- dplyr::bind_rows(results) %>%
mutate(log_sl_x1 = as.factor(round(log_sl_x1, 1)))
# PLOT
ggplot(results, aes(x = preydiv_x1, y = (log_rss))) +
geom_line(size = 1, aes(color = log_sl_x1, group = log_sl_x1)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = log_sl_x1, group = log_sl_x1), alpha = 0.3) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()# Plot
ssf_line_1 +
geom_line(data = results, size = 1, aes(x = preydiv_x1, y = (log_rss),
color = log_sl_x1, group = log_sl_x1)) +
geom_ribbon(data = results, aes(ymin=lwr, ymax=upr, x=preydiv_x1,
fill = log_sl_x1, group = log_sl_x1), alpha = 0.3)# grab stationary probabilities
ps <- momentuHMM::plotStationary(hmm2, plotCI= TRUE, return = TRUE)state1 <- ps$preydiv$'state 1' %>% mutate(state = 1)
state2 <- ps$preydiv$'state 2' %>% mutate(state = 2)
state3 <- ps$preydiv$'state 3' %>% mutate(state = 3)
# bind to one data frame
pdat <- rbind(state1, state2, state3) %>%
mutate(state = as.character(state)) # plot
hmmprobs_plot <- ggplot() +
geom_line(data = pdat, aes(x = cov, y = est, color = state)) +
geom_ribbon(data = pdat, aes(x=cov, y=est, ymax=est+se, ymin=est-se, fill = state),
alpha = 0.4, show.legend = TRUE) +
ylab("Stationary state probabilties") +
xlab("Prey diversity") +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
theme_minimal()
hmmprobs_plot